%用于处理三角网格，完善三角网格的基本信息
function [new_tri, N, dt] = basic_info(ele, nodes, index, varargin)
%index: 1 初始网格 2 加密网格
%varargin: 1 {1}: neighbor, 2{1} 加密次数
if index == 2
    vertices = zeros(size(nodes,1),3);
    for i = 1 : size(nodes,1)
        vertices(i,1:2) = nodes(i,:);
    end
    mesh = surfaceMesh(vertices,ele);
    %numIterations = 1;
    subdivide(mesh,"midpoint-split",varargin{1});
    grid_points = mesh.Vertices;
    new_Faces = mesh.Faces;
    new_nodes = zeros(size(grid_points,1),2);
    for i = 1 : size(grid_points,1)
        new_nodes(i,:) = grid_points(i,1:2);
    end
    ele = new_Faces;
    nodes = new_nodes;
    tight_tri = complete_info(ele, nodes);
end
N = size(ele,1);
min_len = 100;
[b_edge_value, edge_node] = compute_edge_gauss();
for i = 1 : N
    tri(i).nodes = [nodes(ele(i,1),:);nodes(ele(i,2),:);nodes(ele(i,3),:)];
    tri(i).sides = zeros(2,2,3);
    tri(i).sides(:,:,1) = [tri(i).nodes(1,:);tri(i).nodes(2,:)];
    tri(i).sides(:,:,2) = [tri(i).nodes(1,:);tri(i).nodes(3,:)];
    tri(i).sides(:,:,3) = [tri(i).nodes(2,:);tri(i).nodes(3,:)];
    tri(i).circucenter = find_circumcenter(tri(i).nodes);
    tri(i).edgelen = zeros(1,3);
    [tri(i).edgelen(1), tri(i).edgelen(2),tri(i).edgelen(3)] = compute_edge_length(tri(i).nodes);
    min_len = min([min_len,tri(i).edgelen]);
    [tri(i).type, tri(i).which_angle] = judge_type(tri(i).edgelen);
    tri(i).normal = zeros(3,2);
    tri(i).normal = get_normal(tri(i));
    tri(i).alpha = max([abs(sum(tri(i).normal(1,:))),abs(sum(tri(i).normal(2,:))),abs(sum(tri(i).normal(3,:)))]);
    if index == 1
        tri(i).neighbor = varargin{1}(i,:);
    else
        tri(i).neighbor = tight_tri(i).neighbors;
    end
    tri(i).edge_gauss = real_and_ref1(tri(i).nodes,edge_node);
    tri(i).neighbor_rf = zeros(4,2,3);
    tri(i).nei_sides = zeros(1,3);
end
dt = 0.1*min_len/sqrt(2);
[left_b, right_b, up_b, down_b, bd] = category_grid(tri,N,[0,1],[0,1]);
ud_re = find_re_cell(tri,up_b, down_b, 1);
lr_re = find_re_cell(tri,left_b,right_b,2);
%填补每个三角单元的邻居
%先处理左右边界
states_edge = zeros(N,3);
for i = 1 : size(lr_re,1)
    neighbor_rf_nodes = find_rf_nodes_lrbd(tri, lr_re(i,:), [0,1]);
    lc = lr_re(i,1);
    ls = lr_re(i,2);
    rc = lr_re(i,3);
    rs = lr_re(i,4);
    tri(lc).neighbor_rf(:,:,ls) = neighbor_rf_nodes(:,:,1);
    tri(lc).nei_sides(nei_map(ls)) = rs;
    tri(rc).nei_sides(nei_map(rs)) = ls;
    states_edge(lc,ls) = states_edge(lc,ls) + 1;
    tri(rc).neighbor_rf(:,:,rs) = neighbor_rf_nodes(:,:,2);
    states_edge(rc,rs) = states_edge(rc,rs) + 1;
end
%再处理上下边界
for i = 1 : size(ud_re,1)
    neighbor_rf_nodes = find_rf_nodes_udbd(tri, ud_re(i,:), [0,1]);
    uc = ud_re(i,1);
    us = ud_re(i,2);
    dc = ud_re(i,3);
    ds = ud_re(i,4);
    tri(uc).neighbor_rf(:,:,us) = neighbor_rf_nodes(:,:,1);
    tri(dc).neighbor_rf(:,:,ds) = neighbor_rf_nodes(:,:,2);
    tri(uc).nei_sides(nei_map(us)) = ds;
    tri(dc).nei_sides(nei_map(ds)) = us;
    states_edge(uc,us) = states_edge(uc,us) + 1;
    states_edge(dc,ds) = states_edge(dc,ds) + 1;
end
for i = 1 : N
    for l = 1 : 3
        if states_edge(i,l) ~= 0 
            continue;
        end
        [neighbor_rf_nodes, nei, ns] = find_rf_nodes_gen(tri,i,l);
        tri(i).neighbor_rf(:,:,l) = neighbor_rf_nodes(:,:,1);
        tri(nei).neighbor_rf(:,:,ns) = neighbor_rf_nodes(:,:,2);
        tri(i).nei_sides(nei_map(l)) = ns;
        tri(nei).nei_sides(nei_map(ns)) = l;
        states_edge(i,l) = states_edge(i,l) + 1;
        states_edge(nei,ns) = states_edge(nei,ns) + 1;
    end
end
%填补一下左右边界处的邻居
for i = 1 : size(lr_re,1)
    lc = lr_re(i,1);
    ls = lr_re(i,2);
    rc = lr_re(i,3);
    rs = lr_re(i,4);
    tri(lc).neighbor(nei_map(ls)) = rc;
    tri(rc).neighbor(nei_map(rs)) = lc;
end
%填补一下上下边界处的邻居
for i = 1 : size(ud_re,1)
    uc = ud_re(i,1);
    us = ud_re(i,2);
    dc = ud_re(i,3);
    ds = ud_re(i,4);
    tri(uc).neighbor(nei_map(us)) = dc;
    tri(dc).neighbor(nei_map(ds)) = uc;
end
%存一下每条边gauss点相对应在邻居单元处基函数的值
for i = 1 : N
    nei_edge_value = neighbor_edge_value(tri,i);
    tri(i).nedge_value = nei_edge_value;
end
new_tri = tri;
end

function edge_rf_nodes = real_and_ref1(nodes, edge_nodes)
edge_rf_nodes = zeros(size(edge_nodes));
Jacobi = [-nodes(1,1)+nodes(2,1), nodes(3,1) - nodes(1,1);
    -nodes(1,2)+nodes(2,2), nodes(3,2) - nodes(1,2)];
for i = 1 : 3
    for j = 1 : 4
        xi = edge_nodes(j,1,i);
        eta = edge_nodes(j,2,i);
        % 计算形函数
        N1 = 1-xi-eta;
        N2 = xi;
        N3 = eta;

        % 计算实际坐标 (x,y)
        x = N1*nodes(1,1) + N2*nodes(2,1) + N3*nodes(3,1);
        y = N1*nodes(1,2) + N2*nodes(2,2) + N3*nodes(3,2);
        edge_rf_nodes(j,1,i) = x;
        edge_rf_nodes(j,2,i) = y;
        % b = [x- nodes(1,1);y-nodes(1,2)];
        % tmp = Jacobi\b;
        % disp(tmp);
        % disp([xi,eta]);
    end
end
end

function [len12,len13,len23] = compute_edge_length(node)
len12 = compute_distance(node(1,:), node(2,:));
len13 = compute_distance(node(1,:), node(3,:));
len23 = compute_distance(node(2,:), node(3,:));
end

function len = compute_distance(node1, node2)
len_square = ((node1(1) - node2(1))^2) + ((node1(2) - node2(2))^2);
len = sqrt(len_square);
end

function [type, which_angle] = judge_type(len)
% 1 锐角， 2 直角， 3 钝角
% which_angle :1 12最长，2 13最长， 3 23最长
[max_len, which_angle] = max(len);
min_len = min(len);
len_sort = sort(len);
tmp_len = len_sort(2);
judge = min_len^2 + tmp_len^2 - max_len^2;
if judge > 0
    type = 1;
elseif abs(jugde) < 1e-13
    type = 2;
elseif judge < 0 
    type = 3;
end
end

function  normal = get_normal(tri)
normal = zeros(3,2);
index = 1;
for i = 1 : 2
    for j = 2 : 3
        if ( i ~= j)
            x_mid = 0.5*(tri.nodes(i,1) + tri.nodes(j,1));
            y_mid = 0.5*(tri.nodes(i,2) + tri.nodes(j,2));
            normal(index,1) = x_mid - tri.circucenter(1);
            normal(index,2) = y_mid - tri.circucenter(2);
            normal(index,:) = normal(index,:)/norm(normal(index,:));
            index = index + 1;
        end
    end
end
if tri.type == 1
    return;
elseif tri.type == 2
    max_edge_index = tri.which_angle;
    switch max_edge_index
        case 1
            node1 = tri.nodes(1,:);
            node2 = tri.nodes(2,:);
        case 2
            node1 = tri.nodes(1,:);
            node2 = tri.nodes(3,:);
        case 3
            node1 = tri.nodes(2,:);
            node2 = tri.nodes(3,:);
    end
    k = (node2(2) - node1(2)) / (node2(1) - node1(1));
    normal(max_edge_index,:) = [1,k]/sqrt(1 + k^2);
    tmp = zeros(1,3);
    index = 1;
    for i = 1 : 2
        for j = 2 : 3
            if(i~=j)
                tmp(index) = normal(i,:)*normal(j,:)';
                index = index + 1;
            end
        end 
    end
    if min(tmp) > 0 
        normal(max_edge_index,:) = [-1,-k]/sqrt(1 + k^2);
    end
elseif tri.type == 3
    max_edge_index = tri.which_angle;
    normal(max_edge_index,:) = -normal(max_edge_index,:);
end
end

function [left_b, right_b, up_b, down_b,bd] = category_grid(tri,N, x_inter, y_inter)
left_b = zeros(N,2);
right_b = zeros(N,2);
up_b = zeros(N,2);
down_b = zeros(N,2);
bd = zeros(1,N);
index_l = 1;
index_r = 1;
index_u = 1;
index_d = 1;
index = 1;
for i = 1 : N
    if min(tri(i).neighbor) == 0
        bd(index) = i;
        index = index + 1;
        for j = 1 : 3
            if tri(i).neighbor(j) == 0
                switch j
                    case 1
                        wh_bd = judge_bd(tri(i).nodes(2,:),tri(i).nodes(3,:),x_inter,y_inter);
                        side = 3;
                    case 2
                        wh_bd = judge_bd(tri(i).nodes(1,:),tri(i).nodes(3,:),x_inter,y_inter);
                        side = 2;
                    case 3
                        wh_bd = judge_bd(tri(i).nodes(1,:),tri(i).nodes(2,:),x_inter,y_inter);
                        side = 1;
                end
                switch wh_bd
                    case 1
                        down_b(index_d,:) = [i,side];
                        index_d = index_d + 1;
                    case 2
                        up_b(index_u,:) = [i,side];
                        index_u = index_u + 1;
                    case 3
                        left_b(index_l,:) = [i,side];
                        index_l = index_l + 1;
                    case 4
                        right_b(index_r,:) = [i,side];
                        index_r = index_r + 1;
                end
            end
        end
    end
end
left_b = left_b(1:index_l - 1,:);
right_b = right_b(1:index_r - 1,:);
up_b = up_b(1:index_u - 1,:);
down_b = down_b(1:index_d - 1,:);
bd = bd(1:index - 1);
end

function wh_bd = judge_bd(node1,node2,x_inter, y_inter)
%1 down, 2 up, 3 left, 4 right
if (node1(2) == y_inter(1)) && (node2(2) == y_inter(1))
    wh_bd = 1;
elseif (node1(2) == y_inter(2)) && (node2(2) == y_inter(2))
    wh_bd = 2;
elseif (node1(1) == x_inter(1)) && (node2(1) == x_inter(1))
    wh_bd = 3;
elseif (node1(1) == x_inter(2)) && (node2(1) == x_inter(2))
    wh_bd = 4;
else
    error("出错，无法匹配到边界");
end
end
